% Compute the Dominguez and Lobato (2004)'s estimator along with confidence interval
% 
% Arguments:
% h must be a function handle from R^JxTheta into R
% Y must be an n*J matrix of data, which includes both dependent and independent variables.
% X must be an n*K vector of conditioning variables, i.e., instruments.
% hdot must be a function handle referring to the gradient of h w.r.t. θ in row vector.
% p is the dimension of the parameter, θ.
% 
% Examples:
% h = @(y, theta) y(:, 1) - y(:, 2) * theta;
% hdot = @(y, theta) -y(:,2);
% p = 1;

function output = Dominguez_Lobato(Y, X, h, hdot, p, varargin)
%--------------------------------------------------------------------------
% Optional inputs
%--------------------------------------------------------------------------
Tstart = tic();
opt = inputParser;
addParameter(opt, 'lb', -100)  % Lower bound for PSO algorithm
addParameter(opt, 'ub', +100)  % Upper bound for PSO algorithm
addParameter(opt, 'alpha_level', 0.1)
parse(opt, varargin{:});
lb = opt.Results.lb;
ub = opt.Results.ub;
alpha_level = opt.Results.alpha_level;
%--------------------------------------------------------------------------
% Estimation
%--------------------------------------------------------------------------
n = size(X, 1);
K = size(X, 2);
% Define a matrix whose (i, j) element is given by 1{X_i<=X_j}, where <= represents vector ordering
ord = zeros(n, n);
for i = 1:n 
    for j = 1:n 
        ord(i, j) = (sum(X(i, :) <= X(j, :)) == K);
    end
end
ord2 = ord' * ord;

% Define the objective function
obj = @(theta) h(Y, theta)' * ord2 * h(Y, theta);

% Particle Swarm Optimization
options = optimoptions('particleswarm', 'Display', 'off', 'UseParallel', false);
theta_hat = particleswarm(obj, p, lb*ones(1, p), ub*ones(1, p), options);
%--------------------------------------------------------------------------
% Variance estimation
%--------------------------------------------------------------------------
den = hdot(Y, theta_hat)' * ord2 * hdot(Y, theta_hat);
num = hdot(Y, theta_hat)' * ord2 * diag(h(Y, theta_hat) .^2) * ord2 * hdot(Y, theta_hat);
% Variance estimator
Vhat = den \ num / den;
Tend = toc(Tstart);
%--------------------------------------------------------------------------
% Output struct
%--------------------------------------------------------------------------
output.theta_hat = theta_hat;
output.Vhat = Vhat;
output.se = sqrt(diag(Vhat));
output.rtime = Tend;
output.upperCI = theta_hat + norminv(1-alpha_level/2) * output.se;
output.lowerCI = theta_hat - norminv(1-alpha_level/2) * output.se;;
end
